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In this paper, the general procedure to solve the General Relativistic Hydrodynamical(GRH) 
equations with Adaptive-Mesh Refinement (AMR) is presented. In order to achieve, the GRH equa- 
tions are written in the conservation form to exploit their hyperbolic character. The numerical 
solutions of general relativistic hydrodynamic equations are done by High Resolution Shock Cap- 
turing schemes (HRSC), specifically designed to solve non-linear hyperbolic systems of conservation 
laws. These schemes depend on the characteristic information of the system. The Marquina fluxes 
. . . with MUSCL left and right states are used to solve GRH equations. First, different test problems 

^^ I with uniform and AMR grids on the special relativistic hydrodynamics equations are carried out to 

^3 ■ verify the second order convergence of the code in ID, 2D and 3-D. Results from uniform and AMR 

^^ ' grid are compared. It is found that adaptive grid does a better job when the number of resolution is 

Cn ' increased. Second, the general relativistic hydrodynamical equations are tested using two different 

r^ , test problems which are Geodesic flow and Circular motion of particle In order to this, the flux part 

r-{ ■ of GRH equations is coupled with source part using Strang splitting. The coupling of the GRH 

h-^ ' equations is carried out in a treatment which gives second order accurate solutions in space and 

00 ■ *™''- 

,__l ; I. INTRODUCTION 

>; 

^^ , Numerical relativity in astrophysics is becoming an exciting area. The general relativistic observational data 

requires integration of general relativistic hydrodynamical equations with traditional tools which need high compu- 
tational power. There are different numerical tools to solve the relativistic equations, such as general relativistic 
^5 ■ hydrodynamics, magnetohydrodynamics, and nuclear astrophysics. 
T^ I Relativity is necessary for describing astrophysical phenomena involving compact objects. Some of these phenom- 

^^ . ena arc X-ray binary systems, coalescing compact binaries, black hole formations, gamma-ray burst, core collapse 
^j I supernovae, pulsars, microquasars and active galactic nuclei. The general relativistic effects must be considered when 
^H. strong gravitational fields are encountered, such as near compact binaries. The significant gravitational wave signal 
I ' produced by some of these phenomena can also be understood in the framework of the general relativistic theory. 
(2J[), Numerical relativity in astrophysics is becoming more exiting area of research due to the large amount of data 

?• ■ taken from high energy astronomical phenomena. Accordingly, the numerical relativistic astrophysics may help to 
.'^ [ understand some of the observations in high energy and gravitational wave astrophysics, fn order to take advantage 
S^ • from observational data, relativistic equations must be solved to understand phenomena in the strong relativistic 
5h ' regions. The numerical solution of strong field astrophysical phenomena requires integration of the Einstein theory 
with other tools of astrophysics which are relativistic hydrodynamics and nuclear astrophysics. The Einstein theory 
describes the dynamics of the space time with a set of the most complicate equations, ft is highly nonlinear and 
thousands of terms are coupled together in these equations. 

The main issues in General Relativistic Hydrodynamic (GRH) are the numerical modeling of flows with large Lorentz 
factors (strong gravitational regions) and strong shock waves. An accurate description of relativistic flows with strong 
shocks is needed for study of important problems in astrophysics. General relativistic effects caused by the presence 
of strong gravitational fields appear to be connected to extremely fast flows in different astrophysical scenarios, such 
as accreting of compact objects, stellar collapse, and coalescing compact binaries (Banyuls et al., 1997). In the last 
ten years, much effort has been addressed to develop accurate numerical algorithms to solve the equations of general 
relativistic hydrodynamics in extreme conditions. The main conclusion that has emerged is that modern algorithms 
exploiting the hyperbolic and conservative character of the equations are more accurate at describing relativistic flows 
than traditional finite difference upwind techniques with artificial viscosity(Norman et at, 1986). To take advantage 
of the conservation properties of the system, modern algorithms are written in conservative form, in the sense that 
the variation of mean values of the conserved quantities within the numerical cells is given, in the absence of the 
sources, by the fluxes across the cell boundaries. Furthermore, the hyperbolic character of the system of equations 
allows one to treat physical discontinuities appearing in the flow consistently (the shock-capturing property) (Hirsch 
et al., 1992). Over the last few decades, many High Resolution Shock Capturing (HRSC) methods have arisen in 
classical computational fluid dynamics (Hirsch et al., 1992). Nowadays, the numerical simulations employ these HRSC 
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methods. These techniques produce highly accurate numerical approximations in smooth regions of the flow and 
capture the motion of unresolved steep gradients without creating spurious oscillations. 

Our goal is to understand how thing formed in universe. Our strategy is the application of known physics as it is 
dictated by the observations of the objects whose origin and evolution we wish to comprehend, using and creating 
advanced computer techniques. Our tactics are to proceed based on three-dimensional numerical hydrodynamics using 
any source for computing methods wherever feasible, implementing the targeted physical processes in sequence, and 
proceeding to increasingly geometrically complex situations. 

II. FORMULATIONS 
A. Description of Hydrodynamic Equations 

In this section, we give a general technical description of the Hydrodynamic equation as they will be used in our 
code development and the analytic description of different problems. 

The spacetime contains flow in 4-momentum. Each particle carries its 4-momentum vector along its world line. 
The particles or fluid have a 4- velocity u'^(/x = 0, 1, 2, 3) at each point which is the 4- velocity of an observer co-moving 
with the fluid at that point. The unit of 4-velocity is unit of speed of light, c. At each point there are three scalar 
functions to define the characteristic of the fluid. The first of these is the rest mass density p{g/cm?). The other 
scalar functions are the specific internal energy e{ergs/g) and the pressure P in the units of pe which is the energy 
density. The pressure is connected to the other variables by the equation of state P = P{p,e). The quantity ph, 
where h is called the relativistic enthalpy, is the total inertia carrying mass energy. 

The time evolution of the fluid is governed by the General Hydrodynamical equations. The first equation is the 
normalization of the 4-velocity, defined as 

u'^Uf, = -1. (1) 

The GRH equations in references Font et ai, 2000 and Donat et ai, 1998, written in the standard covariant form, 
consist of the local conservation laws of the stress-energy tensor T'*'^ and the matter current density J'^: 

Vm T'^" = 0, v^J^ = 0. (2) 

Greek indices run from to 3, Latin indices from 1 to 3, and units in which the speed of light c = 1 are used. Vm 
stands for the covariant derivative with respect to the 4-metric of the underlying spacetime, g^^. 

In order to quantify the fluid or particles which are moving in their world line, the stress-energy tensor T^'^ is 
defined. The stress-energy tensor exists at each event in spacetime. The stress-energy density is a quantity which 
contains information about energy density, momentum density and stress as measured by any and all observers at 
that event. The stress-energy tensor for an imperfect fluid is defined by Misner et al. as 

2Tja^"' + q^'u'' + q•'u^', (3) 

where the scalar functions ij and ( are the shear and bulk viscosities, respectively. This is the stress-energy tensor for 
a viscous fluid with heat conduction. The scalar 6 is defined as 

= u%. (4) 

It describes the divergence of the fluid world lines. The symmetric, trace-free, and spatial shear tensor a'^'^ is deflned 
by 

ct'^" = ^K^/i^" + W'.^h^'') ~ i0/l^^ (5) 

where h^^ = u-yUf^ + g^^ is the spatial projection tensor. The vector g^ is the energy flux vector. 



Eqs.(Pl-|ll describe the motion of imperfect fluid flow for a fixed metric in complete generality. Entropy for imperfect 
fluid is influenced by viscosity and heat flow. 

Defining the characteristic waves of the general relativistic hydrodynamical equations is not trivial with imperfect 
fluid stress-energy tensor. The viscosity and heat conduction effects are neglected. This defines the perfect fiuid 
stress-energy tensor. This stress-energy tensor is used to derive the hydrodynamical equations. Using perfect fluid 
stress-energy tensor, we can solve some problems which are solved by the Newtonian hydrodynamics with viscosity, 
such as those involving angular momentum transport and shock waves on an accretion disk, etc. Entropy for perfect 
fluid is conserved along the fluid lines. The stress energy tensor for a perfect fluid is given as 

T^"' = phuf'u'' + Pgt"". (6) 

A perfect fluid is a fluid that moves through spacetime with a 4- velocity u^ which may vary from event to event. It 
exhibits a density of mass p and isotropic pressure P in the rest frame of each fluid element, h is the speciflc enthalpy, 
defined as 

h = ! + €+-. (7) 

P 

Here e is the speciflc internal energy. The equation of state might have the functional form P = P{p^ e). The perfect 
gas equation of state, 

P={T- l)pe, (8) 

is such a functional form. 

The conservation laws in the form given in Eq.(|2Jl are not suitable for the use in advanced numerical schemes. 
In order to carry out numerical hydrodynamic evolutions such as those reported in Font et ai, 2000, and to use 
HRSC methods, the hydrodynamic equations after the 3-fl split must be written as a hyperbolic system of flrst 
order flux conservative equations. The Eq.© is written in terms of coordinate derivatives, using the coordinates 
(x^ = t,x^ ,x^,x^). Eq.® is projected onto the basis {n'', (^7)'^}, where n^ is a unit timelike vector normal to a 
given hypersurface. After a straightforward calculation, we get (see ref.Font et ai, 2000), 

dtU + d,F'^S, (9) 

where dt = d/dt and di = d/dx"^ . This basic step serves to identify the set of unknowns, the vector of conserved 
quantities [/, and their corresponding fluxes FiU). With the equations in conservation form, almost every high 
resolution method devised to solve hyperbolic systems of conservation laws can be extended to GRH. 

The evolved state vector U consists of the conservative variables (D, S'j, r) which are conserved variables for density, 
momentum and energy respectively; in terms of the prinutivc variables {p,v^,e), this becomes(Font et ai, 2000) 

. fD\ ( ^Wp \ 

U^\ SA^\ ^phWS . (10) 

V ^ / V V^iphW ^P-Wp) ) 

Here 7 is the determinant of the 3-mctric 7^ , Vj is the fluid 3- velocity, and W is the Lorentz factor, 

W = av° = {1- -i^jv'v^)-^'^. (11) 

The flux vectors F' are given by (Font et ai, 2000) 

F' = I a{{v^ - ^(3^)S, + 77^,5]} | . (12) 



The spatial components of the 4-velocity m* are related to the 3-velocity by the following formula: u' = W{v'^ — P^ jo). 
a and /3' are the lapse function and the shift vector of the spacetime respectively. The source vector S is given by 
(Font et al, 2000) 






(13) 



where F" is the 4-dimensional Christoffel symbol 



^% == -^9°''^ {d 1.9^13 + 9^5p/3 - 3/35^^). 



(14) 



B. Spectral Decomposition and Characteristic Fields 

The use of HRSC schemes requires the spectral decomposition of the Jacobian matrix of the system, dF^/dU. The 
spectral decomposition of the Jacobian matrices of the GRH equations with a general equation of state was reported 
in Font et ai, 2000. They display the spectral decomposition, valid for a generic spatial metric, in the x-direction, 
{pF^ /dU)\ permutation of the indices yields the other two directions. 

They started the solution by considering an equation of state in which the pressure P is a function of p and e, 
P = P{p, e). The relativistic speed of sound in the fluid Cg is given by(Font et ai, 2000) 



C 



dP 

dE 



X 
h 



Pk 



(15) 



where x = dP/dp\e, k = dP/de\p, S is the entropy per particle, and E = p + pe is the total rest energy density. 
A complete set of the right eigenvectors [t^;] and corresponding eigenvalues A^ along the x-direction obeys 
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\f!i\ = Xi [fl] , i = 1,...,5. 



The solution contains triply degenerate eigenvalues (Font et ai, 2000), 



(16) 
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A^ = A2 = A3 = aiJ - p . 
The other eigenvalues are given by (Font et ai, 2000) 
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± 



A set of linearly independent right eigenvectors spanning this space is given by (Font et ai, 2000) 
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hWiK-pC^) 



,Vx,Vy,Vz,l 



hW{K-pC^) 



(19) 



WVy, h{-1^y + 2W'^Va,Vy), h{-fyy + 2W'^VyVy), 



h{-ty, + 2W^VyV,), VyW{2Wh " l) 



(20) 



WV^hij,, + 2W^V,V,), h{jy, + 2W^VyV,), 



H-f,, + 2W^v,v,),v,W{2Wh - 1) 



(21) 



and 



hhWiv, 



w^ - (A± + /3^)/a 



,hWvy, hWvz, 

-[T 
-1 



(22) 



Here, the superscript T denotes the transpose. 

Now, the eigenvahies and eigenvectors are computed in the y and z-directions using the information in the x- 
direction. Let M^(f/) be the matrix of the right eigenvectors, that is, the matrix having as columns the right 
eigenvectors with the standard ordcring(Banyuls et al., 1997), M^ = (r_,ri, r2, ^3, r+). To obtain the spectral 
decomposition in the other spatial directions q = 2, 3(= y, z), it is enough to take into account the following symmetry 
relations: 

1) The eigenvalues are easily obtained by carrying out in equations (|17|l and H18|l a simple substitution of 
indices x hy q. 

2) The matrix M-^ in the x-dircction is permuted according to M.'^ = Pqi(M.^), where the two operators Pqi 
act on their arguments by the following sequential operations: 

• to permute the second and (q + l)th rows, 

• to interchange indices 1 = a; with q. 

From Ibancz et al, 1998, the left eigenvectors in the a;-direction are : 



lo,i 



W 



h-W 
Wvy 

-w 



(23) 



h^- 



-JzzVy + JyzVz 
v'^ilzzVy -lyzVz) 

lzz{^ - VxV') + JxzVzV'' 

-7^2(1 - VxV'') - JxzVyV'' 

-IzzVy + lyzVz 



(24) 



i0,3 — 



h£_^ 



-lyyVz + IzyVy 

V''{,lyyV^--1,yVy) 

-lzy{^ - VxV"^) - ^xyVzV"" 

7yy{l - V^V"") + J^yVyV"^ 

-lyyVz + IzyVy 



(25) 



(±1) 



A^ 



r,,(l - KAX) + (2/C - \)vi{w^v^e - ^xccvn 



r,j,(i - icAi) + (2/C - i)vi{w^vye - r^y«") 

r,,(l - ICAl) + (2/C - l)VUW^v'e - ^xzv^) 
(1 - /C)h7w^ + V^(M^2^^ _ r^^)] - icw^vw 
The variables used for the left eigenvectors in the x-dircction are defined as follows: 
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Al = ~X±+(3'' , A = A/a , /3^ = /3^/c 



/C = 



K — C% 



, k = k/ p 



Cl = v,-Vi , Vi 






A% = 



ryXX _ yXj^X_ 



1-Ai^ v^v^ , Ai-A^^^ w"(c| - c;) 



(c| - c|) + (i;v^ - ^^v;) = 



A"" = h-'WilC - 1){CI - Cl)^'' 



r = r,, - ^v^v^ 



^ — dCt^ij — ^xx^ XX ~r '~ixy^ xy i ^xz^ xz 



r.. 


— lyylzz - 


"fyzlzy 


-L xy 


= lyzlzx ~ 


~ ^yx^zz 


r.z 


= lyxlzy - 


~ lyylzx 



71) 

Now, the left eigenvectors of the Jacobian matrix are computed in the y and z-directions using the left eigenvectors 
in the x-direction. Let L^([/) be the matrix of left eigenvectors, that is, the matrix having as rows the left eigenvectors 
with the standard ordering. The eigenvalues of the left eigenvectors are the same with the eigenvalues of the right 
eigenvectors. To obtain the spectral decomposition in the other spatial directions 9 = 2, 3(= y, z), it is enough to take 
into account the following symmetry relations: 

The matrix Ad^ in the a;-direction is permuted according to L** = VqiiJ-^), where the two operators Vqi act on their 
arguments by the following sequential operations: 

• to permute the second and [q + l)th columns, 

• to interchange indices 1 = a; with q. 

III. NUMERICAL METHOD 
A. Marquina Fluxes 

Approximate Riemann solver failures and their respective corrections (usually adding a artificial dissipation) have 
been studied in the literature(Quirk, 1994). Motivated by the search for a robust and accurate approximate Riemann 
solver that avoids these common failures, Shu et al., 1989 have proposed a numerical flux formula for scalar equations. 
Marquina flux is generalization of flux formula in Ref.Shu et al., 1989. In the scalar case and for characteristic 
wave speeds which do not change sign at the given numerical interface, Marquina's flux formula is identical to 
Roe's flux(Toro, 1999). Otherwise, scheme is more viscous, entropy satisfying local Lax-Friedrichs scheme(Shu et al., 
1989). The combination of Roe and Lax-Friedrichs schemes is carried out in each characteristic field after the local 
linearization and decoupling of the system of equations. However, contrary to other schemes, the Marquina's method 
is not based on any averaged intermediate state. 

The Marquina fluxes(Font et al., 2000) with MUSCL left and right states are used to solve the 3-D relativistic 
hydro equation. In Marquina's scheme there are no Riemann solutions involved (exact or approximate) and there are 
no artificial intermediate states constructed at each cell interface. 

To compute the Marquina fluxes, first, the sided local characteristic variables and fluxes are computed . For the left 
and right sides, the characteristic variables are 



LP{Ui)-Ui, wP^LP{Ur)-Ur (27) 



and the characteristic fluxes are 



^P^LP{Ui)-FiUi), ^P = LP{Ur)-F{Ur). (28) 

where the number of conservative variables p = 1..5. Ui and Ur are conservative variables at the left and right sides, 
respectively. U'{Ui) and LP{Ur) are the left eigenvectors of the Jacobian matrice, dF^/dU. 

The left and right fluxes are defined depending on the velocities of the fluid for each specific grid zone. The 
prescription given in Ref.Donat et al., 1998 is as follows. 

For all conserved variables p ~ 1, ..m 

if Xp{U) does not change sign in (if ( \p{Ul) x Xp{Uji) > 0)), then 
if \p{Ui) > then 



else 

$P = $P 
end if 
else 

ap = max[/er(c/,,(7,) |Ap([/)| 
$P = 0.5($f + akwf) 
$P =0.5($P + ap<) 
end if 
where Ap is an eigenvalue of the Jacobian matrix and, 

afe=max{|Ap([/0|,|Ap(t/,)|}. (29) 

The numerical flux that corresponds to the cell interface separating the states Ui and Ur is then given by Ref.Donat 
et ai, 1998: 

F''{Ui, Ur) = ^($^rP([/0 + <^P_rP{Ur)). (30) 

p=i 

Marquina's scheme can be interpreted as a characteristic-based scheme that avoids the use of an averaged interme- 
diate state to perform the transformation to the local characteristic fields. 

In carrying out Marquina's scheme, we have to compute intermediate states and the Jacobian matrix of the states 
at each cell interface. So we need to know the left and right states, Ul and Ur, at each interface. To construct the 
second-order scheme, the MUSCL left and right states are used, 

uf = lii(O) = Uj - -Aj 

uf = u,iAx) = u, + -A,. (31) 

To avoid the appearance of oscillations around discontinuities in MUSCL-type schemes, slope limitcrs in the recon- 
struction stagc(Hirsch et ai, 1992, Toro, 1999) are used .A common TVD limiter is based on the minmod function 
(Hirsch et ai, 1992). The standard minmod slope provides the desired second-order accuracy for smooth solutions, 
while still satisfying the TVD property. Minmod slope function is given as 

A,; = minmod(Ui - Ui_i, Uj+i - Uj), (32) 

where the minmod function of two arguments is defined by: 

a if \a\ < \b\ and ab > 
minmod(a, b) — -^ b if |6| < \a\ and ab > (33) 

if a6 < 0. 

B. Strang Splitting 

HRSC methods are applicable without source term. We now turn our attention to handling the source and flux 
terms. To solve the general relativistic hydrodynamical equations numerically Strang splitting is used, in which Eq.Q 
is split into two parts. The first is given by the flux tcrms(Bona et al., 1998) 

dU dF"" dpy dF'- _ 
dt dx dy dz 

The evolution operator of Ea. (|34|l can be described by an F, so the solution of Ea. l34() becomes 



U(t + At) = FU(t) (35) 

The second has the source terms 

^=S. (36) 

With the evolution operator for Ea. l|36|l . S, the solution of Ea. l|36|) is 

U{t + At)^SU{t) (37) 

In the numerical solution, this splitting is performed by a combination of both flux and source operators. Denoting by 
U{At) the numerical evolution operator for the system in a single time step, the evolution operator for solution(Toro, 
1999) is written as: 

Uit + At) =. Si^)FiAt)S{^)Uit). (38) 

The resulting update of U{t + Ai) is second-order accurate in space and time if and only if both the operators S and 

F areO(Aj;2,At2). 

This kind of splitting allows easy implementation of the different numerical treatments of the principal part of the 

system without worrying about the sources of the equations. It also allows choosing the best scheme for each type of 

problem. 

C. Updating the Equations in Time 

The HRSC methods will be used to solve the hydrodynamical equations. This requires that the hydrodynamical 
equations are written in a flux form. 

Solution of the relativistic hydrodynamical equations using MUSCL left and right states with Marquina fluxes gives 
second-order accurate solutions in space. In order to get solutions that are also second-order accurate in time, the 
solution is evolved in the following way. 

To solve Ea. (|38|l first, the solution of Eq.(j37|) should be computed, which is called a system of Ordinary Differential 
Equations (ODEs). Using the two stage Runge-Kutta Tnethod{ToTO, 1999), which gives a second-order accurate 
solution in time, 

Ki = AtS{f\U"-) 
A'2 = AtS{t" + At, U" + R'l) 

[/"+i = [/« + l[Xi+A'2]. (39) 

Explicit methods are much simpler to use than implicit methods. The latter require the solution of non-linear 
algebraic equations at each time step and are therefore much more expensive. If ODEs which have no spatial 
variations, F(U) = 0, are solved by the implicit method, then there will be no stability restriction on the time step. 
Since an explicit method is used to solve the ODEs here, the time step is limited by the Courant condition (Hirsch et 
ai, 1992), 

At = C „ , ^^, , (40) 

(|w"| + C,)„,ax ^ ' 

where C is called Courant number. It can be < C < 1. Ea. (|40|l means that fastest wave at a given time travels for 
at most one cell length Aa; in the sought time step At. 

Second, the solution of Ea. H34|) is computed using a predictor-corrector time update method, as follows. 
Initially, the solution at the half time step is computed. 
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2A2/ 
At 



\\fji,j+^,k — {f"')i,]-^,k) 



ark.ku (/")«.,,fc-i)- (41) 



Then, the sohition using the full time step is computed and new fluxes {f"^^), which arc functions of C/, 



»,i,fe ' 



At 

'i,3,k — ^ij.fc ~ 'Xx^^'' ' >i+^^3^k - \J ' )i-^d 

At 



urn = vi^M - ^((r+^).+i.,,fe - (r+^),_i ,fe) 



|J((r+^k,>+i -(/"+* W-i)- (42) 

In order to solve general relativistic hydrodynamical equation in 3 — 13, the following procedure is followed. First 
Ea. (|36|l is solved using second-order Runge-Kutta method, Ea. (|39|l . at half time step. Later Eg . (|34|l is solved using 
predictor-corrector method, Eas. (|41l42|l . at the full time step. Finally, the source term is solved another half time step 
using the Rungc-Kutta method. At the end of the calculation, a second-order accurate solution for GRH equation is 
obtained in space and time. 

IV. ADAPTIVE-MESH REFINEMENT 

Adaptive mesh refinement is needed to solve some problems in numerical relativity, such as coalescing compact 
binaries, accreting compact objects, etc. A key difficulty in calculating accurate waveforms from these simulations 
lies in the fact that the length and time scales for these problems can vary by an order of magnitude or more within 
the computational domain. Adequate resolution is needed to model strong field regions. Adaptive techniques are thus 
crucial for success in modeling such systems. 

A. Paramesh Package 

The Paramesh package (MacNeice et ai, 2000) are used to make to implement AMR and parallelization in the 
code. The package runs on any parallel platform which supports either the shmem(shared memory) or MPI libraries. 
This package is a set of FORTRAN 90 routines designed to enable Adaptive Mesh Refinement (AMR) for a parallel 
distributed memory machine. 

The basic AMR strategy is that the computational domain is covered with a hierarchy of numerical sub-grids. 
These sub-grids form the nodes of a tree data-structure and are distributed among the processors. When more spatial 
resolution is required at some location, the highest resolution sub-grid covering that point spawns child sub-grids (2 in 
ID, 4 in 2D and 8 in 3D), which together cover the line, area or volume of their parent, but now with twice its spatial 
resolution. They refer to the block with the highest refinement level at a given physical location as a 'leaf block'. All 
sub-grids have identical logical structure (i.e., the same number of grid points in each dimension, the same aspect 
ratios, the same number of guard cells, etc ). They differ from each other in their physical sizes and locations, and 
in their list of neighbors, parents and children. The package assumes that the application will use logically Cartesian 
grids (i.e., grids that can be indexed i, j, k). It works for ID, 2D and 3D problems. 

The package provides routines which perform all the communication required between sub-grids and between 
processors. It manages the refinement and de-rcfincment processes, and can balance the workload by reordering the 
distribution of blocks among the processors. It also provides routines to enforce conservation laws at the interfaces 
between grid blocks at different refinement levels. Flux conservation at different refinement level interfaces will be 
explained in section llVBI . 

The top panel of Fig^ shows a 6 x 4 grid structure in 2D on each block. The numbers assigned to each block 
represent the block location in the three structure seen in the bottom panel of Fig^ Notice that in Paramesh the 
refinement level is not allowed to jump by more than 1 level at any location of the spatial domain. 
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B. Flux Conservation 

Conservation of conserved variables in the hydrodynamical code is an important physical result. To guarantee 
conservation of variables on the AMR grid, some special operation at block boundaries of different refinement levels 
are needed. Conservative hydrodynamics codes require that the fluxes entering or leaving a grid cell through a common 
cell face shared with the 4 refined neighbors, equal the sum of the fiuxes across the appropriate faces of the 4 smaller 
cells, Fig|21 (MacNeice et ai, 2000). This is achieved by taking the average of the fine gird fiuxes and copying them 
into the coarse grid to correct fiuxes at different refinement levels. For more information about this fiux conservation, 
check Ref.MacNeice et at, 2000. 

C. Refinement Criteria 

A refinement criteria needs to be defined to analyze any problem using Paramesh package. We have chosen to use 
a refinement criterion based on the first derivative of density, which has different dynamics during the evolution for 
the shock tube problem. When density is refined or dercfined, the other variables, such as pressure and velocity, will 
be refined or dercfined simultaneously. Two different functions are used to define refinement criteria, which are given 
as follows: 



Tiix,y) 



f \dx dy 



(43) 



and 



T2{x,y) 



(/(*,J)-/C» + 1,J + 1)) 



(44) 



where Eas. H43|) and (|44|l are the first derivative and diagonal derivative of a function /, respectively. The function / 
can be chosen any variable in the simulation depending on the change of dynamics during the numerical simulation. 
Refinement and derefinemcnt of blocks depend on refinement criteria and parameters. If the value of the refinement 
functions, Ea. (|43|l or H44|) . are bigger than refinement parameter, ctore specified by user based on problem, in any zone 
in the same block, that block will be refined. If the value of the refinement functions are smaller than derefinemcnt 
parameter, ctode, that block will be dercfined. Refinement parameter has to be bigger than derefinemcnt parameter. 
These refinement and derefinemcnt criteria are checked on each timestep during the numerical simulation. The 
refinement criteria and parameters play an important role in AMR. Note that both refinement and derefinemcnt in 
Paramesh are restricted by the requirement that the grid can change by only one refinement level at each block. 



RESULTS FROM SPECIAL RELATIVISTIC SHOCK TUBE PROBLEM 



The Riemann shock tube problem (Hirsch et ai, 1992, Font et ai, 2000) is used as a test case. In this problem, 
the fiuid is initially in two different thermodynamical states on either side of a membrane. The membrane is then 
removed. In the test problem, it is assumed that the fluid initially has pl > PR, where the subscripts L and R refer 
to the left and right sides of the membrane. With this assumption, a rarefaction wave travels to the left, and a shock 
wave and contact discontinuity travel to the right. The fluid in the intermediate state moves at a mildly relativistic 
speed which is w = 0.725c to the right. Flow particles accumulate in a dense shell behind the shock wave compressing 
the fluid and heating it up to values of the internal energy much larger than the rest-mass energy. Hence, the fluid is 
extremely relativistic from a thermodynamical point of view, but only mildly relativistic dynamically. The Riemann 
shock tube is a useful test problem because it has an exact time-dependent solution, and tests the ability of the code 
to evolve both smooth and discontinuous flows. In the case considered here, the velocities are special relativistic and 
the method of flnding the exact solution differs somewhat from the standard non-relativistic shock tube. 

The initial values of the fluid variables for ID test are; pi = 10.0, pr = 1.0, ui = 0.0, u^ = 0.0, pi ~ 13.3,pr = 
0.66.10"^. The computational domain is < a; < 1 and at f = the membrane is placed at x = 0.5. We use the 
perfect fluid equation of state Eq.© with F = 5/3. 
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TABLE I: Comparing our code Li norm errors with the other hydrodynamic code results from 3D. 



LI norm error for shock tube problem 


Code 


npts 


ii(p) 


Li{p) 


Li{v) 


OUR 
CODE 


32'^ 
64^ 
128^ 


0.188 
9.813E-2 
5.439E-2 


0.1798 
9.199E-2 
4.595E-2 


1.82E-2 
6.178E-3 
3.86E-3 


Font, Miller 
Suen, Tobias 


128^ 


9.23E-2 


7.98E-2 


9.66E-3 


GENESIS 
Aloy et ai, 2000 


40^ 
60^ 
80^ 
100^ 
150=* 


l.lE-1 
9.8E-2 
9.2E-2 
7.0E-2 
4.8E-2 


8.0E-2 
5.2E-2 
4.5E-2 
3.7E-2 
2.5E-2 


0.9E-2 
l.lE-2 
l.lE-2 
7.0E-3 
5.0E-3 



TABLE IL Convergence rate from 3-D special relativistic hydro code. 



Convergence rate for the shock tube problem 


npts — 2 * npts 


r{p) 


r{p) 


r{v) 


32-^ - 64'^ 
64^ - 128^ 


0.94 
0.851 


0.9668 
1.001 


1.563 
0.68 



Uniform Solution of 3D Shock Tube 



The 31? special relativistic shock tube problem is solved using the initial data which is given above. The MUSCL-left 
and right states and the minmod slope function with the Marquina fluxes are used. This method gives a second-order 
accurate solution for the 31? special relativistic hydrodynaniical equations. The Li norm errors and convergence 
rates are computed and our results are also compared with the existing literature. In order to compare the SD 
numerical with ID analytic solution for the shock tube, the initial discontinuity is set up along the diagonal of the 
computational domain. The 128'^ zones are used with interval of length l/VS in each directions. The diagonal of the 
cube therefore has unit length. The maximum time, t = 0.4, is used to compute analytic and numerical calculations 
and the numerical data is extracted along the main diagonal of the cube. 

Table ^ shows the Li-norm errors of different hydrodynamical quantities. These are the density, pressure and 
velocity. We present our results, those from Font et at, 2000, and GENESIS code(Aloy et ai, 2000) in Table |l| all 
of these are 3D simulations. Our code and the code of Ref.Font et ai, 2000 use the Marquina method with MUSCL 
left and right states which give a second-order accurate solution. The Genesis code uses Marquina method with 
the PPM reconstruction procedure which gives third-order accurate solutions. Table |l] shows that the results are 
comparable. Another difference between the Genesis and our code is that we compute Li-norm errors at t = 0.4 
but they compute at i = 0.5. Marquina method with PPM left and right states is computationally more expensive 
than our code. Because it has more reconstruction process and it needs to more evolution step to get the third-order 
accuracy in time. The PPM interpolation algorithm described in Colella et ai, 1984 gives monotonic conservative 
parabolic profiles of variables within a numerical zones. In the relativistic version of PPM, the original interpolation 
algorithm is applied to zone averaged values of the primitive variables (p, v,p), which are obtained from zone averaged 
values of the conserved variables U. 

The convergence factors are shown in Table ^ It should approach 1 when the higher order methods are used. 
Fig. O shows the result obtained for different numbers of points with this method. The solid line represents the exact 
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solution while different symbol represents the numerical results for different resolution for the density. As can be seen 
from Fig. |31that when the number of resolutions are increased, numerical solution converge to analytic one. 

B. AMR Results 

We have first carried out simulations of the 2D shock tube using 3 levels of refinements. To do this problem, 
initial discontinuity has been placed along the diagonal axis of the computational domain (a square). The evolution 
of the system along the direction orthogonal to this diagonal is purely one dimensional problem. It allows us to 
compare one-dimensional analytic solution with the numerical results. The results of shock tube problem with AMR 
are compared with analytic solution at each grid zone of the computational domain. The error in each grid zone and 
the refinement levels are plotted in Fig^ The finest level is replaced around r ~ 0.8 where shock wave is produced. 
Eg 1431 is used as a refinement function. Refinement parameters for this problem are: ctore(refinement parameter) = 
0.035, ctode(derefinement parameter) = 0.005. 

The refinement function in Eq^|is used to solve the same problem with the same refinement parameters of Fig0] 
The numerical results are given at Fig|5l The results show that Fig|5lis more derefined regions and use less number 
of resolutions. 

In this part of test problem same refinement function is used in Fig[3but refinement parameters for this problem 
are: ctorc(refinement parameter) = 0.035, ctode(derefinement parameter) = 0.007. 

It is seen from Figs 0]|S1 and that refinement levels are different for each cases. FigElhas more derefine regions than 
the others. If the more derefine regions are observed, this indicates that less number of points is used during the 
evolutions. The highest refinement level is replaced at discontinuity. As a conclusion, AMR grid structure depends 
on refinement functions and parameters which need to be defined carefully for different problems. 

Here, the Riemann shock tube problem is solved in 2D using seven different refinement levels with the best re- 
finement parameters which are ctore (refinement parameter) = 0.035 and ctode(derefinement parameter) = 0.007. In 
Figd density is displayed on the x — y plane at the time when the shock tube structure is well-developed. The AMR 
run of the shock tube problem creates a nice grid structure that follows the density gradients in the problem well. 
In Fig|Hl the error function for each grid zone is given with refinement level to see the behavior of the adaptive grid 
during the evolution. We have also looked at the movie for this problem to see the behavior of refine ad derefine grid 
during the calculations. The finest grid follows the shock wave and some derefine regions are produced at contact 
discontinuity which is a constant state. 

It is interesting to compare the computing time for AMR and uniform runs in 31?. Fig|51 illustrates the time per 
processor using different resolution for AMR and uniform runs. The times for the uniform and AMR cases are the 
same for low resolutions, 16'^ and 32^, because the whole computational domain is refined to the maximum level at 
initial time. It means that there is no AMR in these cases. But the computing times for the resolutions 64'^ and 
128'^ in AMR are less than for uniform grids. The computing time for lower resolution is affected by overhead in the 
AMR routines. Saving the computational cost in this test problem for 128^ is approximately 30%. Increasing the 
size of computational domain, the number of refinement levels and the number of resolutions increase saving of the 
computational cost in the shock tube problem. Increasing these futures allow us to see more lower levels grids and 
more derefine region when shock wave moves to left and rarefaction wave move to right. 

Fig. ^1 depicts the fine grid volume fraction in the AMR simulation with 128"^ zones. Initially, less than 15% of 
the computational volume was covered by the fine grid. The fine grid volume fraction increases as the shock wave 
propagates to the left and the rarefaction wave propagates to the right. Although dcrefinement in the computational 
domain occurs around the origin at t = 0.25 the total number of grid cells continues to increase as the shock and 
rarefaction waves move to the right and left, respectively, as shown in Fig. ^1 The computational volume is more 
derefined than refined around t = 0.3. The numbers of refined and derefined zones depend on the refinement criteria 
and parameters used in numerical simulations. 

In these kinds of test simulations, the saving in computational cost of the AMR run over a unigrid run are modest 
since the wavelength of the wave, A is comparable to the dimension L of the computational domain. In the simulation 
of the real astrophysical problem, it is expected that wavelength is much smaller than the dimension of computational 
domain, resulting in significant computational savings. For a full simulation of binary coalescence or accretion disk 
around the compact objects, one will need to handel multiple scales to evolve dynamics of the sources and propagation 
of the waves that develop. While the sources are expected to require high resolution, particularly near regions of shocks 
or shear layers, the waves can likely be handled with lower resolution. 
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VI. GENERAL RELATIVISTIC HYDRODYNAMICAL TEST PROBLEM 



111 this section, we first introduce the Schwarzschild geometry in spherical coordinates to define sources in the general 
relativistic hydrodynamical equations. Then the accretion disk problems are solved , which have been analytically 
analyzed, to test the full GRH code in 2D on the equatorial plane. 

A. Schw^arzschild Black Hole 

The Schwarzschild solution determined by the mass M gives the geometry in outside of a spherical star or black 
hole. The Schwarzschild spacetime metric in spherical coordinates is 



ds^ 



- 1 )dt^ +1 )-^dr^ + 

r r 

r'^de'^ + r'^sin' 



(45) 



It behaves badly near r = 2M; there the first term becomes zero and the second term becomes infinite in Eg . H45(l . 
That radius r = 2M is called the Schwarzschild radius or the Schwarzschild horizon. 
The spacetime metric for this line element is as follows: 



5m^ 






V 



(1 












^)-^ 











r2 











2 ■ 2 



(46) 



The lapse function and shift vector for this metric is given bellow: 



jj"^ = 0.0, /?" = 0.0, /3* = 0.0, 
a=(l-^)V^. 



(47) 



B. The Source Terms For Schw^arzschild Coordinates 

The gravitational sources for the GRH equations are given by Eg. ljl^l) . In order to compute the sources in 
Schwarzschild coordinates for different conserved variables, Ea. H13|l can be rewritten as. 



S 



\ 



( 



(48) 



It is seen in Ea. (|13|) that the source for conserved density, Z?, is zero but the other sources depend on the components 
of the stress energy tensor, Christoffel symbols, and 4-metric. After doing some straightforward calculations, the 
sources can be rewritten in Schwarzschild coordinates for each conserved variable with the following form: 
The source for the momentum equation in the radial direction is 



1 



-^r^7^<^^{T"drau+T^''drgr 



The source for the momentum equation in the 9 direction is 



(49) 
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«V7r^''ff..r^, = \a^T**deg^^. (50) 

The source for the momentum equation in the (j) direction is 

a^T^^'-guaTl^ = 0.0. (51) 



The source for the energy equation is 



a^iT'^'dra ~ aT^'g^'drgu). (52) 

The non-zero components of the stress-energy tensor in Schwarzschild coordinates can be computed by Eq.®; they 
are 





T" = ph 


w^ 

^+^5" 
a^ 


rprr 


=: phW^il 


ff + Pgrr 


rj^ee 


= /9/iVK2(i 


fif + PgO' 


rp<t><f> _ 


= phW^iv 


'I'Y J^Pg'l"i> 




rptr 


= ph V . 

a 



(53) 



C. Geodesies FIoavs 

As a general relativistic test problem, an accretion of dust particles onto a black hole is solved. The exact solution 
for pressureless dust is given in Ref.Hawley ei.aZ, 1984. This problem is evolved in 2D in spherical coordinates at 
constant 6 = 7r/2, which is the equatorial plane. Accordingly, the spatial numerical domain is the (r, (p) plane. The 
2.4 < r < 20 and < < 27r are used for the computational domain. The initial conditions for all variables are 
chosen to have negligible values except the outer boundary (r = 20Af ) where gas is continuously injected radially 
with the analytic density and velocity. They are formulated as: 



and 



2^^ L 2M . ^. 

r V r 



(55) 



where W is the Lorentz factor and is given by 

here M is mass of black hole 

Throughout the calculation, whenever values at the outer boundary are needed, the analytic values are used. The 
code is run until a steady state solution is reached using outflow boundary conditions at r = 2.4, inflow boundary 
conditions at r = 20 and the periodic boundary in the direction are used. No dependence on (p is found throughout 
the entire simulation. 
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TABLE III: Li norm error and convergence factors are given for different resolutions. 



Convergence Test 


^ of points 


# of time step 


I/i norm error 


Convergence factor 


32 


1 


1.675E - 5 




64 


2 


4.0811-B - 6 


4.105 


128 


4 


9.6456£' - 7 


4.23 


256 


8 


2.1372E' - 7 


4.51 



Fig. ^] depicts the rest-mass density p, absolute velocity v = (v'^Vi)^^^ and radial velocity v^ as a function of 
radial coordinate at a fixed angular position. The numerical solution agrees well with the analytic solution. Some 
convergence tests with the SR Hydro code are conducted. They confirmed that it is second order convergent. The 
convergence test are also made on the this problem to test the behavior of the source terms in the GR Hydro code. 
The analytic values of the accretion problem are compared with numerical results. The convergence results are given 
at Tabl dllll It is noticed that code gives roughly second order convergence. 

The conservation form of the general relativistic hydrodynamical equations arc solved. FromEl the general form 
of a conservation law is expressed by stating that the variation per unit time of this conserved equation within the 
volume is zero. Then it is expected that conserved variables must be conserved to machine accuracy, ~ 10~^^, at each 
time step. The results of conservation variables in the numerical test problem show that these variables are conserved 
to machine accuracy. 



D. Circular Motion of Test Particles 

The circular motion of a fluid on the equatorial plane will be simulated. To do this, the circular flow is set up with 
negligible pressure in the equatorial plane, in which angular velocity at each radial direction r is given by: 



(1-^) 



(57) 



It is called a Keplerian velocity. 

The analytic expectation of circular motion is that, the last stable circular orbit for a particle moving around a 
Schwarzschild black hole in a Keplerian disk is at r = QM {M is mass of black hole) . Therefore, the gas or particles 
fall into the black hole, if their radial position is less than 6M. When their radial position is bigger than 6M, they 
should rotate in a circular orbit. Here, this problem is simulated and the numerical solution is compared to the 
analytic expectations. It tests the code with sources in the (j) direction. 

In order to simulate this problem, the computational domain is chosen to be 3A/ < r < 20M and < < 27r. 
The computational domain is filled with constant density and pressureless gas, rotating in circular orbits with the 
Keplerian velocity and zero radial velocity. In Fig. ^1 the radial velocities of the gas vs. radial coordinate at different 
times are plotted. It is numerically observed that the gas inside the last stable orbit, r = 6Af , falls into the black hole, 
while gas outside the last stable orbit follows circular motion with the Keplerian velocity as it is expected analytically. 
Fig. ^1 shows the density of the fiuid vs. radial coordinate at different times to see the behavior of the disk. It is 
seen that gas falls into the black hole for r < 6M. Accordingly, the numerical results from our code are consistent 
with the analytic expectations. 
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VII. CONCLUSION 

111 this paper the Eulerian type numerical sohition of general relativistic hydrodynamic equations constructed for 
relativistic astrophysics is presented. This code is constructed for general spacetime metric with lapse function and 
shift vector. It is capable of solving special relativistic problem and general relativistic problem coupling with the 
Einstein and hydrodynamic equations. This paper explains how general relativistic hydrodynamic equation can be 
solved with AMR. 

In the first part of the paper, the GRH equations and the numerical method are given. GRH equations are written 
in the conservative form and the numerical method and Paramesh package is applied to them with some appropriate 
way to get working code with AMR. 

Then, we have presented and tested solution of three-dimensional relativistic hydrodynamic equations using MUSCL 
scheme with Marquina fluxes. The results are obtained for problems involving relativistic flows, and strong shocks. 
These problems are solved using uniform and adaptive grid. Adaptive-Mesh Refinement is getting one of important 
technique in numerical relativity to solve the real astrophysical problems, such as accretion disks around the compact 
objects, and coalescing of the compact binaries. Adaptive-Mesh simulations pose many complex design problems, and 
a variety of different techniques that exist to solve these problems. The AMR is used to solve relativistic equations 
and to present some results from our numerical calculations. We solved the shock tube problem using different 
refinement levels to see the behavior of adaptive grid and also to compare AMR and uniform results with analytic 
solution. Results are converging to analytic solution when the number of resolutions are increased. Computing time 
used to solve AMR and uniform grid to get same accuracy are also compared. Results show that when the number 
of resolutions are increased AMR uses much less computing time than uniform grid. 

Finally, The full general relativistic hydrodynamical equation are tested using 2 — D test problems which are called 
free falling gas onto black hole and the circular motion of particles around the black hole. The numerical results from 
these test problems agree with analytic expectation and converge to analytic solution when the number of resolutions 
are increased. 

As a conclusion, AMR has permitted us to address problems of astrophysical interest which are really hard to 
solve with uniform grid. Because these problems are dominated by local physics, this code has been designed to 
handle highly relativistic flows. Hence, it is well suited for 2D and 3D simulations of astrophysical problems. First 
astrophysical results will be presented in a forthcoming paper for two-dimensional accretion disk problem around the 
black hole. 
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FIG. 1: Representation of adaptive-mesh refinement grid and tree structure. Figure is taken from MacNeice et at, 2000. 
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FIG. 2: Flux conservation is fixed at different refinement level boundaries. Figure is taken from MacNeice et al., 2000. 



21 



CO 

I 0.5 







o^iH|iHo l iiti l o | o>$»B | o l ^ l QMcp « fc^^ 



Shock Tube 



t = 0.4 




analytic 
32x32x32 
64x64x64 
128x128x128 



jj | ; l $B{!^:p306|]tjJli 






0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

r 



FIG. 3: The analytic and numerical solutions of the relativistic shock tube problem for 3-D. 
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FIG. 4: Plotting the error and the refinement levels in each zone Eg. 1431 is used as the refinement function with parameters: 
ctore (refinement parameter) = 0.035, ctode(derefinement parameter) = 0.005 
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FIG. 5: Plotting the error and refinement levels in each zone. Eg. 1441 is used as the refinement function with parameters: 
ctore (refinement parameter) = 0.035, ctode(derefinement parameter) — 0.005 
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FIG. 6: Plotting the error and refinement levels in each zone. Eg. 14411 is used as the refinement function with parameters: 
ctore (refinement parameter) = 0.035, ctode(derefinement parameter) = 0.007 
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FIG. 7: Density is plotted in the x-y plane. Each square represents the block structure and each block has 8 points in a 
each direction. Eg. 1441 is used as the refinement function with refinement parameters: ctore(refinement parameter) = 0.035, 
ctode(derefinement parameter) — 0.007. 
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FIG. 8: We plot the error and refinement levels in each zone to see the behavior of the refinement levels with errors. Eg. 14411 
is used as the refinement function with refinement parameters: ctore(refinement parameter) = 0.035, ctode(derefinement 
parameter) — 0.007. 
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Total numbers of processor used for all run are 64 
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FIG. 9: The computing time for AMR and uniform grid runs of the 3-D shock tube problem. The red hne is the time per 
processor for AMR and the black line is the time per processor for uniform grid runs. 
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FIG. 10: (a): The fraction of the volume covered by the fines grid for AMR run of the 3D shock tube problem, (b); the density 
of the shock tube problem at different times. 
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FIG. 11: The analytic solutions (red solid lines) with the numerical solutions (black circles) using 256 zones in the radial 
direction for thermodynamical variables, p, v, and v^ . 
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FIG. 12: Radial velocity vs. r is plotted. The radial velocity of fluid in the disk stays zero for r > 6M, during the evolution. 
r = 6M is called the last stable circular orbit in a Keplerian disk. 
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FIG. 13: Density vs. r is plotted. Density of the fluid is plotted at different times using different colors. The matter falls into 
the black hole while r is less than 6M. 



